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ABSTRACT 

Oh- 

Context. The light curves of variable stars are commonly described using simple trigonometric models, that make use 
of the assumption that the model parameters are constant in time. This assumption, however, is often violated, and 
■ consequently, time series models with components that vary slowly in time are of great interest. 

Aims. In this paper we introduce a class of data analysis and visualization methods which can be applied in many 
different contexts of variable star research, for example spotted stars, variables showing the Blazhko effect, and the 
spin-down of rapid rotators. The methods proposed are of explorative type, and can be of significant aid when performing 
a more thorough data analysis and interpretation with a more conventional method. 

Methods. Our methods are based on a straightforward decomposition of the input time series into a fast "clocking" 
periodicity and smooth modulating curves. The fast frequency, referred to as the carrier frequency, can be obtained 
. from earlier observations (for instance in the case of photometric data the period can be obtained from independently 

^1^' measured radial velocities), postulated using some simple physical principles (Keplerian rotation laws in accretion disks), 

or estimated from the data as a certain mean frequency. The smooth modulating curves are described by trigonometric 
polynomials or splines. The data approximation procedures are based on standard computational packages implementing 
simple or constrained least-squares fit -type algorithms. 
^ ■ Results. Using both artificially generated data sets and observed data, we demonstrate the utility of the proposed 

methods. Our interest is mainly focused on cases where multiperiodicity, trends or abrupt changes take place in the 
variable star light curves. 

Conclusions. The presented examples show that the proposed methods significantly enrich the traditional toolbox for 
^ ■ variable star researchers. Applications of the methods to solve problems of astrophysical interest will be presented in 

' the next papers of the series. 

' Key words, stars: variables: general - methods: data analysis 

1. Introduction of a global model, a certain set of local variability elements 

111 r- 1 1 11 1 rr postulatcd, and then their change in time is followed. 

^ The light curves of variable stars can display very different probably the best known method of this kind is the wavelet 

modes of behaviour. They can be quite regular, often even transform 

>■ periodic, or relatively fast changes can be detected. The ^ ' , , i . , , . r- 

methods to deal with such observed data can roughly be ^ ^o^ continuously observed signals the Fourier transform 



O 



X 



, divided into two classes. The first class encompasses tools "based methods and different types of time-frequency trans 

H that are based on Fourier analysis methods, assuming that ^^^^ are a well-developed part of the traditional data pro- 

the light curve can be decomposed into a set of periodic ^^^^^'^8 work flow (see e.g. Poularikas {ed.),mQ)- The sit- 

(harmonic) components. The components themselves can ^^tion is different in the case of astronomical data The ir- 

then be interpreted as manifestations of, e.g., different pul- ^^g^lf'^ observing patterns and long time gaps m the mea- 

sational modes (as in asteroseismology), rotation, orbiting, ^^^^^d series significantly complicate the analysis (see e.g. 

and so forth. The local behaviour of such multicomponent review by Schwarzenberg-Czerny, 2003, and references 

curves can be quite complex. For instance, it is quite easy t^^erem or Chapter 5 in Aerts et al.20ia). For the first type 

to build simple examples where destructive interference be- methods, the computed frequency spectra often show 

tween two nearby components with nearly equal frequen- spurious peaks due to the regularities in the gap structure 

cies allows some third component to dominate the observed f ^he input data. The second class of methods, on the other 

curve. The detection of such local regularities often lead to ^^^d, can only be used for the short, more or less continu- 

the usage of the second type of analysis tools, usually re- subsets of data. 

ferred to as the time-frequency analysis methods. Instead In this paper we propose a relatively new method, at 

least in the astronomical context, to deal with irregularly 

Send offprint requests to: J. Pelt spaced data. We start from the simple physical assump- 
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tion that even if an object is observed to exhibit irregu- 
lar or semi-regular variability, some regular process can be 
identified in the background. The best example of this is 
stellar rotation, normally quite stable over long periods of 
time, the semi-regularity of the observed light curve be- 
ing caused by inhomogeneities, such as dark spots, on the 
stellar surface. Another relevant setting is provided by ac- 
cretion disks, where the Keplerian rotation of the accretion 
elements is more regular than the observed outfiow of radi- 
ation. Using this assumption, we decompose the observed 
light curves into two components: a rapidly changing car- 
rier tracing the regular part of the signal, and its slowly 
changing modulation. In this way, we construct a method 
combining both aspects of the two classes of methods pre- 
sented above. 

The proposed method is certainly not a new method to 
seek hidden periods from data with gaps and irregularities. 
Instead, we try to build continuous models for irregularly 
observed data using a priori known frequency. Proper visu- 
alization of the model can then help us to interpret the data 
at hand and to reveal possible physical effects. The method 
is of explorative nature and can be used as a starting point 
before building more sophisticated statistical, numerical or 
dynamical models. 

The layout of the paper is somewhat non-traditional. 
Instead of commencing with a lengthy introduction and 
comparison of the new proposed method to various time 
series analysis methods that have been used in variable star 
research, we start with introducing the basic elements of the 
analysis method in Sect. [2l where we also give a detailed 
description of the algorithms used. Sect. [3] is dedicated to 
demonstrate the suitability of the proposed method in var- 
ious astronomical settings, using both artificially generated 
and observational data sets. This allows the reader to get 
appreciation of the new method at work. In Sect. |4] we will 
give an overview of the related developments, discuss the 
new method in the light of previous work, and in Sect. [5] 
discuss the domains of its applicability. The full-scale ap- 
plication of the method to longer sets of observational data 
will be presented in the forthcoming papers in the series. 



2. Carrier fit method 

The notion of the carrier frequency is well known from com- 
munication theory. Very often this is just a single harmonic 
tone with a fixed frequency vq- The carrier is modulated 
by a signal of interest and transmitted to receiver. In the 
receiver the modulated signal is demodulated and useful 
information is recovered. The simplest model for such a 
process is just a waveform 



f{t) = a(t) cos(27rti^o) + ^(f) sin(27rti^o), 



(1) 



where the carrier wave is modulated by two smooth low- 
frequency signal components a{t) and h{t). In some cases 
the signal can be represented in an analytic complex form 



Ut) = A{t)e^^i^^'\ 



(2) 



which allows the explicit tracking of the time dependent 
amplitude A{t)^ the phase (f{t)^ and most importantly, the 
instantaneous frequency 



d(p{t) 
dt 



The real astronomical signals very often contain also higher 
harmonics of the basic frequency. Due to this we introduce, 
from the very beginning, a more complex model of the form 



K 



f{t) = ao{t)^Yl cos(27rt/ci/o) + bk{t) sm{27Ttkiyo)) .{4) 

k=i 

By including the term ao(t), we allow the mean level of the 
signal to be different from zero. We also assume that the 
total number of harmonics, is not very large. 

2.1. The algorithm 

To describe the details of the carrier fit method we start 
from the simplest model of a variable star light curve: 



fit) = a{t) cos(^) + b{t) sin(^), 



(5) 



where Pq is the period (measured in days) corresponding 
to the carrier frequency vq = I/Pq (measured in cycles 
per day), and the two functions a{t) and b{t) describe slow 
changes of the curve. Here we define a slow process as one 
for which the characteristic changes of the functions a and 
b occur on time-scales significantly longer than the period 
of the carrier oscillation, Pq. 

To build approximation curves for the observed data 
we firstly need to have a proper value for the carrier pe- 
riod Pq, the determination of which is discussed in detail in 
Sect. 12.51 Secondly, certain analytical or numerical models 
for the both slow modulating functions are required. These 
models themselves depend on certain sets of parameters, 
as discussed in Sects. \T2\ and 12.31 for the chosen two repre- 
sentations. After reliably determining Pq and formulating a 
suitable model for the modulators, the remaining task is to 
find the model that best fits the observations, and retrieve 
the corresponding parameters. This data fitting process it- 
self can be formulated as a standard least-square approxi- 
mation procedure. This procedure can, in general, involve 
both linear and non-linear parameters. Because of the over- 
all complexity of non-linear fitting, we usually restrict our 
analysis to the linear regime. 

Below we demonstrate the carrier fit tool using two dif- 
ferent methods for smooth approximation of the modula- 
tors: trigonometric polynomials and splines. There are cer- 
tainly many other possibilities - essentially every known 
method of data smoothing from literature can also be used 
in this context. 



2.2. Trigonometric approximation 

We start from the relatively transparent method of trigono- 
metric approximation. This allows us to explicitly control 
the frequency domain behaviour of our algorithm. 

Let the time interval [tmin^tmax] be the full span of 
our input data. Then we can introduce a certain period 
D = C X {tmax — tmin) for which the coverage factor C 
is larger than unity (typically C = 1.1 — 1.5). Using the 
corresponding frequency, = 1/P^, we can now build a 
trigonometric (truncated) series of the type: 



(3) 



a{t) = Co + ^ [ci cos{2iTtlvD) + sin{2iTtlvD)), 



(6) 
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and 



bit) 



(7) 



^=1 



to be used as models for smooth modulating curves. The 
coverage factor C must be larger than the data span to 
avoid artefacts due to the built-in periodicity of our trigono- 
metric models. Too large coverage factor leads to redundant 
parametrization. Very often the building of the trigonomet- 
ric model starts from the rounding of the full span period, 
and checking that the coverage factor remains in reason- 
able limits. According to our definition of a slow process, 
the period D must be significantly longer than the carrier 
period Pq. 

Using basic trigonometric identities it is now easy to 
show that the full expansion of the carrier fit model will 



contain cosine and sine terms of the form 



Ii^d)) where / = 0, . . . , L. Proper expansion coefficient esti- 
mates must be computed for every term in the series for the 
fixed carrier frequency and "data frequency" ud] this is 
a standard linear estimation procedure and can be imple- 
mented using standard mathematical (statistical) packages. 

In the beginning of our discussion about the algorithm, 
we started off with the simplest light curve model, Eq. (j5j), 
including only the lowest harmonic /c = 1. To generalise to 
the full case, Eq. (|4]), including also higher harmonics, we 
will construct trigonometric representations ak{t) and bk{t) 
for each overtone of the carrier kuo^ k = 1^ . . . ^ K . The time 
dependence of ao{t) is also modelled using special trigono- 
metric polynomial. If all polynomials use the same number 
of harmonics L and we approximate separate cycles by in- 
harmonic model, then the overall count of linear parameters 
to be fitted is TV = (2 x L + 1) * (2 x + 1). 

The actual choice of the representative parameters K 
and L depends on the particular object we are working 
with. The number of tones, depends on the complexity 
of the phase curves. The choice of L is constrained by the 
longest gaps in the time series. In principle, it is also pos- 
sible to compute certain formal constraints using standard 
statistical regression theory (see Draper & Smith 1998). 
Very often the inspection of phase plots which are com- 
puted with different parameter sets can also be useful. 

2.3. Spline approximation 

The trigonometric models of the modulation curves are the 
most useful when the actual behaviour of the modulations 
is relatively continuous, i.e. without sharp transients. The 
spline approximations are more flexible from this point of 
view. 

To build spline approximations we need to divide the 
full data time span into, say L, intervals. For every interval 
with index / = 1, . . . , L, and for each overtone of the carrier 
/c = 1, . . . , in, we define two local sets of cubic curves 



m=0 

and 

3 

bki{t) = ^ bkimt^, 

m=0 



(8) 



(9) 



to be used as local modulators for the different tones of 
the cosine and sine components of the carrier wave. The 
complete analytical model for the data will then consist of 
cosine components 



dkimt'^ cos(27rHi^o), 
and sine components 



(10) 
(11) 



There will altogether be 2 x AKL such components. To 
take into account the general level of the signal, we need 
to include also the unit component, after which the total 
number of modes to be fitted into data will be 8KL + 1. 
The splines are not just a sample of local polynomials, 
but must obey certain restrictions. In our simple case, we 
demand that the approximation curves themselves, their 
first and second derivatives should be continuous. As a re- 
sult we need to complement our set of fitting equations 
with 2 X 3K{L — 1) linear constraints. The total number 
of free parameters will then he N = 2KL + + 1. To 
perform the constrained least-squares fitting we used rou- 
tines from the AS 164 -package written by Stirling (I 198ip . 
The implemented algorithms are based on the so-called 
Givens rotations (see Gentleman I1973p and allow to solve 
weighted least-squares equations. The constraints are effec- 
tively taken into account by assigning infinite weights to 
the constraining equations. 

2.4. Construction of an analytic signal 

One of the useful properties of the carrier fit method is that 
it allows easily to compute Hilbert transforms for the model 
waveforms. Using the Bedrosian theorem (Bedrosian ll963p . 
we can easily compute for term 

fit) = a{t) cos(27rti/o) + b{t) sin(27rti/o), 
its Hilbert transform 



h{t) = a{t) sin(27rti^o) — ^(^) cos(27rti^o)- 



(12) 



By definition, the modulators a{t) and b{t) are smooth 
functions, or in spectral terms, their spectra are concen- 
trated around the zero frequency. The carrier frequency is 
significantly higher, and therefore the modulators can be 
removed from the transform. The complex analytic signal 
(Gabor ri946p can now be constructed as 



w{t) = f{t)+jh{t), 

its amplitude as function of time 



Ait) = ^p(t) + hHt), 

and the instantaneous frequency 



h'f - f'h 
MP + h^) 



(13) 



(14) 



(15) 



can be computed. For the details we refer the reader to the 
paper by Vakman (.1996) . 

Unfortunately the analytic signal approach is useful 
only for waveforms which are either completely or nearly 
harmonic. For more general curves the instantaneous values 
for amplitude, phase and frequency can be formally com- 
puted, but their interpretation is not straightforward. 
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2.5. Estimation of the carrier frequency 

The raost important input parameter for the carrier fit 
method is naturally the carrier frequency vq, or the corre- 
sponding period Pq = l/^o itself. Very often a useful value 
is already known from previous data analysis steps. Here 
are some typical cases: 

— Estimate from previous literature is at hand. 

— Orbiting period is known from eclipses. 

— Rotation period is computed from radial velocity data 
(say, in fully synchronized binaries). 

— The period is estimated by averaging its seasonal values 
(Blazhko effect stars, spotted stars). 

— The central frequency of the band in the Fourier spec- 
trum can be used (Quasi-Periodic Objects). 

In most cases, therefore, an initial value can be easily de- 
termined. The first trial can be improved later if the phase 
plots show significant upward or downward trends. 

In the case when we do not have a proper value for car- 
rier frequency we need to estimate it from the data. Next 
we formulate a method to compare different carrier frequen- 
cies, and describe how a rough initial value can be improved 
to eliminate any trends in phase plots. 

The most natural choice to measure the quality of the 
carrier frequency value is to inspect the stability of the re- 
sulting light curve - with the optimal carrier frequency, the 
most stable behaviour should be obtained. For this pur- 
pose we can use the interpolated carrier fit trial waveform. 
Even if built using a very rough preliminary estimate for 
the frequency vq , it approximates the input waveform in its 
entirety and is a continuous function. Let us fix a certain 
trial frequency ut somewhere around the preliminary car- 
rier frequency vq . We can divide the interpolated waveform 
into intervals of the length Pt = and construct local 

phase curves 

y^M = f{ti + PT<p), (16) 

where ti = to -\- iPt-, i = 0, 1, ... , Nt — 1 is a sequence of 
starting points of the different intervals of the interpolated 
curve /(t), phase (j) runs from to 1, and Nt is the total 
number of the intervals for a particular trial period. The 
phase curves change if we move along time (later we will 
plot such time-dependent curves as phase diagrams). We 
can now measure the rate of change of the phases along 
time using a statistic 

^ i=o i 

Loosely speaking, we estimate the mean square of the first 
derivative of the phase diagram along time. The optimal 
value for carrier frequency should minimize our statistics 
(see the discussion in Sect. 13.91 and Fig. [22]) . 

After computing the value for an optimal carrier fre- 
quency, we can recompute the curve estimates and repeat 
the optimal frequency search, and thereby obtain a refined 
value for the optimal carrier frequency. In actual compu- 
tations such a procedure takes only a couple of iterations 
to converge. Using this approach, uncertainties in the ini- 
tial value of the carrier frequency do not matter, as the 
method ultimately finds the optimal and consistent carrier. 



However, one problem can be identified with this proce- 
dure. The optimal carrier frequency, as defined above, is 
essentially only a representative device. Sometimes it may 
be hard to find a physical interpretation to to the result- 
ing carrier value computed from the data. Only if we can 
assume that behind the seemingly erratic variation of oscil- 
lation parameters there is a certain mean flow which can be 
characterized by a certain mean frequency, we can accept 
the solution as physically meaningful. Fortunately, in most 
practical cases, the carrier frequency is available from the 
very beginning or physical considerations lead to a proper 
interpretation easily. 



2.6. Visualisation 

We use the following scheme to visualize our results. First 
we calculate a continuous curve estimate, /(t), from the 
randomly spaced and gapped data set, using a carrier-based 
least-squares fitting scheme. This approximation is contin- 
uous and does not contain gaps, and therefore allows us to 
get a smooth picture of the long-term behaviour. Below we 
will discuss why and when we can trust such a gap-filling 
scheme. Next we divide this continuous curve into segments 
with a length of the carrier period Pq = Throughout 
the paper, we measure periods in the units of days, 
and frequencies in the units of cycles per day. We 
then normalize each segment so that the approximating val- 
ues span the standard range of [—1, 1]. After normalization, 
we stack the segments along the time axis. To enhance the 
obtained plot, we somewhat extend every segment along 
phases, so that the actual display is wider (along phases) 
than a single period. The normalization is a relevant part 
of our procedure because it helps to grasp phase informa- 
tion we are interested in (trends, drifts, flip-flops etc). In 
the lower part of the plots we often depict the distribu- 
tion of the actual observational time moments in the form 
of "bar code". The black stripes correspond to the peri- 
ods during which at least one observation is available. The 
bright stripes correspond to gaps. This method of visuali- 
sation allows to verify that the model fits into the data and 
not into the gaps. If the number of nodes or harmonics L 
for the modulation model curves is chosen properly then 
the phase plots do not reveal underlying timing structure. 



2. 7. Local time-frequency analysis 

The obtained continuous, 'de-gapped' curve f{t) obtained 
by the carrier fit method can be used as an input to different 
additional analysis procedures. Now that a fully sampled 
data set has been obtained, the standard approaches are 
applicable, due to which many different choices are possible. 
For instance, we can compute local Fourier spectra for slid- 
ing time windows, or compute various wavelet transforms 
(see Poularikas (ed.), 2010. for additional possibilities). The 
corresponding time-frequency plots can then be combined 
with phase diagrams to get full grasp of the changes in the 
underlying processes involved. 

In particular, we use stacks of normalized, local in time, 
power spectra to reveal trends in frequency structure of 
artificial or observed curves. 
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Fig. 1. V-band photometry of LQ Hya, JD 2447881 - 
2452053. The longest time gap in the data set is 284 days. 
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Fig. 2. V band photometry of MW Lyr, JD 2453920 - 
2454052. 



3. Examples and applications 

The carrier fit method is quite straightforward in its for- 
mulation and uses rather standard tools for its implemen- 
tation. However, its full utility can be appreciated only by 
looking at some actual applications. To build up the re- 
quired intuition, we first illustrate the method applying it 
to very simple analytical light curve models which are sam- 
pled at time points taken from real observations, and after 
that, to real astronomical data sets. In the current context, 
however, we totally ignore the physical background of the 
effects observed in each example star - our goal here is to 
show the versatility of the method. 

3.1. Timing sequences and data sets 

To get realistic long-term photometric timing sequences we 
use a compilation of the V-band photometry for the vari- 
able star LQ Hya (Berdyugina et al. i2002;) . The used sub- 
set of 1627 points (displayed in Fig. [T]) covers the inter- 
val JD 2447881-2452053 and contains rather long gaps (up 
to 284 days). In the most of the examples below we use 
only time points from this data set and calculate artificial 
time-dependent process values using analytical formulae. 
To model observational errors we add to each analytical 
time series value a noise component e{t)^ which we gener- 
ate as a normally distributed random variable with stan- 



Fig.3. V-band photometry of FK Com, JD 2451144 - 
2451380. This is the only data fragment where the flip-flop 
phenomenon is observed continuously. 

dard deviation of 5% from the overall series amplitude. This 
level of scatter is quite characteristic for reasonably good 
photometry, corresponding roughly to the signal-to-noise 
ratio of 20. 

To demonstrate the suitability of the new methods to 
analyse the so-called Blazhko effect, we use the RR Lyr 
-type variable star MW Lyr. The data and detailed photo- 
metric solutions for the extensive photometry can be found 
in Jurcsik et al. (2008). As an input for reanalysis we used 
a well-sampled subset of the observations spanning over the 
interval JD 2453920 -2454052, containing 3216 points (see 
Fig.©. 

The usability of the method in the context of the flip- 
flop phenomenon (e.g. Jetsu et al. I1993P is demonstrated 
using a subset of the extensive photometry of the rapidly 
rotating variable star FK Com, described in Korhonen et 
al. (2007) and depicted in Fig. [3l The data set contains 
968 points and spans over the time interval JD 2451144 - 
2451380. 

The periods used as analytical model parameters are 
also taken to match some values that can be derived from 
observations; we return to their actual deduction from the 
data in the next papers of the series. 

3.2. Mismatched carrier frequency 

Very often the appropriate carrier frequency uq is obtained 
from information which is available before the actual data 
processing. Being it rotation, orbiting or previously deter- 
mined frequency, we need to confirm or re-evaluate it. To 
do this, a model with smoothly varying coefficients is fitted 
into the data, using the procedure described in Sect. 12.51 
If the actual frequency is stable, but differs from the used 
carrier frequency we will see it from the time-dependent 
phase diagram. In the following example we use artificially 
generated data, but use the time-points from actual obser- 
vations of LQ Hya described above. For instance, if the 
input model is just a single harmonic with a frequency 
ly = 0.149009 = 1/6.711001, magnitude offset and an added 
noise term £{t) (see Sect. 13. ip 

f{t) = 12.0 + cos(27ri/t) + e{t), (18) 

then using a carrier value = 0.148714 = 1/6.724333, we 
will get a time-dependent phase diagram depicted in Fig. S) 
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Fig. 4. Time-dependent phase diagram, i.e. the hght curve 
amphtude profile over phase (y-axis) plotted as function of 
time (x-axis), for the mismatched carrier frequency exam- 
ple presented in Sect. 13.21 In this plot, a slightly too low 
carrier frequency of = 0.148714 is used, with K = 1 
model and L = 8 modulation harmonics. Time points are 
obtained from real V-band photometric observations of LQ 
Hya; the bar-code in the bottom of the plot indicates when 
data has been available (black) and the gaps (bright). For 
visualisation details see Sect, f 
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Fig. 5. Time-dependent local spectra computed for the in- 
terpolated curve /(t), to recover the correct carrier fre- 
quency uq. See Sect. [372l 



To stress that the actual time points used in the analysis 
are taken from real observations with long gaps, we often 
plot a stripe-like illustration of the observing times on the 
lower part of the phase diagrams, as explained in detail in 
Sect. 12.61 With the refined-search method for the carrier 
frequency, we can recover the actual period in the input 
data set, after which the tilted stripes in the figure become 
strictly horizontal. 

To analyse time-dependent frequency spectra, we use 
the same continuous estimate, /(t), but now we use it as an 
input to a standard time-frequency softwar^B- In our par- 
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Fig. 6. Time-dependent local spectra computed for the 
original data points without interpolation, to be compared 
with Fig. [5] 



ticular case we compute frequency spectra for sliding time 
windows. After computing each local in time spectrum we 
again normalize it to span the standard range [0, 1]. In this 
case the normalization helps us to track minor shifts in fre- 
quencies even if the amplitudes of the separate components 
change drastically. In this way our visualization method en- 
hances the aspects of the variability which we are looking 
for. The result of this analysis is shown in Fig. [5l 

In principle it is possible to compute local frequency 
spectra for the original (not interpolated) data set. In Fig. [6] 
we depict the results of this kind of transform. ^From this 
plot it is evident how the gaps in the data tend to dis- 
tort the overall image. The changes in the spectra can be 
(mistakenly) counted as real modulations. 

3.3. Abrupt change in the period 

The next interesting feature in variable star light curves is 
a possible jump in the frequency. Applying the same time 
points from observations as in the previous example, such 
a process can be described with a simple time series model 



fit) 



_ j cos(27rt/Pi) ■ 
cos(27rt/P2) ■ 



-e{t),t < 3000, 
-e{t),t > 3000, 



(19) 



^ ISDA - Irregula rly Spaced Data Analysis, available at 
http://www.aai. ee/^ pelt /soft. htm 



where Pi=6. 747017, P2=6. 701800, e{t) is noise term de- 
scribed above (Sect. 13. ip We choose the carrier frequency, 
i^o, as a mean of the two model frequencies, and apply the 
carrier fit method to build a continuous curve /(t); the 
results are depicted in Fig. [71 As a comparison, local time- 
frequency analysis is made, the results of which are plotted 
in Fig. [HI The smooth appearance of the two plots is a re- 
sult of the proper approximation scheme. The gaps of the 
input data are no longer visible, as they are properly filled 
in. This facilitates the visibility of the continuous patterns 
in the changing phase and period. However, if to compare 
the two plots, we see quite a relevant difference: the tran- 
sition time from one period to the other is much longer 
for the time-dependent spectra if compared with the phase 
plot. To explain the difference, we need to take into ac- 
count the fact that the time resolution of the transients 
differs for the two analysis methods. The smoothness of the 
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Fig. 7. The effect of an abrupt change in the period 
(Sect. 13. 3p in the time-dependent phase diagram. The pe- 
riod before the jump is Pi = 6.747017, and after it P2 = 
6.7018004. The used carrier is Pq = 6.724333, the number 
of model harmonics, K = 1, and the number of modulation 
harmonics, L = 8. 




Time (days) 

Fig. 8. The effect of an abrupt change in the period in the 
time-dependent local spectra; to be compared with Fig. [71 

phase plots is determined by the smoothness of the modu- 
lating coefficients a{t) and and can be adjusted using 
different fitting schemes. For time-dependent spectra the 
smoothness is controlled by the width of the sliding trans- 
form window. To obtain sufficient resolution in frequency 
domain we lose resolution in time domain. This is a well- 
known and ubiquitous uncertainty relation, and is common 
for all time- frequency analysis methods. 

3.4. Two beating waves 

Next we consider an example, where an interference be- 
tween two nearby frequencies occurs. Again, the same time 
points with gaps as in the previous examples are used, and 
an artificial time series model of the form 

f{t) = cos(27rt/Pi) + cos(27rt/P2) + e{t), (20) 

where Pi=6. 747017, P2=6. 701800, and e{t) is noise term 
(see Sect. 13. ip . is analysed using the carrier fit method to 
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Fig. 9. The effect of two beating periods. Pi = 6.747017 
and P2 = 6.701800, in the time-dependent phase diagram. 
The carrier used is Pq = 6.724333, the number of model 
harmonics K = 1, and the number of modulation harmon- 
ics P = 8. See Sect. [3711 

build a phase diagram, plotted in Fig. [9l where a clear 
chessboard pattern can be seen. This can be understood 
by considering the trigonometric identity 

A-\- B A — B 
cos(A) + cos(P) = 2 cos( — - — ) cos( — - — ). (21) 

The interference between the two nearby frequencies can be 
looked upon as a certain high-frequency main tone modu- 
lated by a low-frequency beat tone. In our case the beat pe- 
riod is ^ 2000 days. Because we are dealing here with 100% 
modulation, the sign of the beat wave changes after every 
1000 days, and the overall image of the time-dependent 
phase looks like a chessboard. These abrupt changes in the 
phase are somewhat unexpected if we take into account that 
the model contains only two simple harmonics. The carrier 
fit method and the used visualization technique allows us 
to reveal this natural flip-flop in the clearest way. 

The continuous approximation curve f{t) can be fur- 
ther analysed by using the local time-frequency analysis 
method. The results of such an analysis with two different 
sliding windows are shown in Figs. [10] and [Til The first of 
these two figures, in which a sliding window of the length 
600 days was used, demonstrates a very serious problem 
of the time- frequency analysis. From the form of the input 
model, we expect that two parallel lines along time are re- 
covered. Instead, we see that the plot is dominated by a sin- 
gle maximum at a period which corresponds to the mean of 
the two model frequencies, and only around the locations, 
where the modulating beat wave changes sign, the spec- 
trum splits into two separate frequencies. Even then the 
local maxima are not placed at the correct periods. Such 
a weird picture is due to the usage of a relatively short 
sliding time window (600 days) if to compare with beat pe- 
riod (2000 days). Even if we compute the local spectra for a 
four times longer time window (as in Fig. [11]), there is still a 
slow wave seen in the positions of the two frequencies. This 
example illustrates a typical error made in time series anal- 
ysis, where the original data is divided into segments, and 
each segment is analysed separately, after which conclusions 
are made according to the obtained local periods. Were the 
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Fig. 10. The effect of two beating periods in the time- 
dependent local spectra, sliding window length is 600 days; 
to be compared with Fig. [9l 
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Fig. 11. The effect of two beating periods in the time- 
dependent local spectra, sliding window length is 2400 days; 
to be compared with Figs. [QlandfTOl 
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Fig. 12. The effect of a linear trend in frequency, as de- 
scribed in Sect. 13. 5[ in the time-dependent phase dia- 
gram. The used model parameters for the carrier fit read 
= 0.148714, = 1, and L = 8. 
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Fig. 13. The effect of a linear trend in frequency in the 
time-dependent local spectra. To be compared with Fig.[T2l 



local spectra computed without a preceding continuous in- 
terpolation, the results could be even more misleading. In 
this particular case we consider the time-dependent phase 
plot a more appropriate analysis and visualization method. 

3.5. Linear trend in the frequency 

As our next example of the carrier fit method, we define ^ 
the input data as follows 

f{t) = cos(27rt((l-A(t)x0.148214+A(t)x0.149214))+£:(t), (22) 



S 0.1490 
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where X{t) = ^ ^ ^rnin moves linearly from to 1 during 
the time span of the input data [tmin^tmax] and e{t) is a 
noise term (see Sect. 13. ip . 

As seen from Fig. [121 the phase diagram for a curve 
with a changing frequency is rather complex, and it is not 
easy to interpret the actual form of the variability in the 
input data. In this case, however, the local time- frequency 
analysis reveals quite explicitly what is going on, as shown 
in Fig. [HI 



Fig. 14. The instantaneous frequency estimated using the 
Hilbert transform scheme of the carrier fit utility (thin line) 
overplotted with the analytical model (thick line) for the 
example presented in Sect. 13.51 Parameters for the carrier 
fit model: uq = 0.148714, K = 1 and L = 8. 
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Fig. 15. The effect of the interference of non-harmonic 
waves with ui = 0.148214 and U2 = 0.149214 in the 
time-dependent phase diagram for a carrier fit model us- 
ing uo = 0.148714. See Sect. IsTgI 

For this model, we can also apply an instantaneous fre- 
quency estimation scheme based on the Hilbert transform 
to construct an analytic signal, as described in Sect. 12.41 In 
Fig. El the computed run of the instantaneous frequency 
is displayed. Due to the noise component and gaps in the 
data set, the recovery is not exact. Also, the boundary ef- 
fects in the beginning and end of the data are clearly visi- 
ble, resulting from the extra freedom for the corresponding 
parts of the fitting curve. The results for this example also 
demonstrate an interesting aspect of time-frequency analy- 
sis. ^From Eq. (|22|) we can see that for the last time point, 
the frequency in the input data is 0.149214. The instanta- 
neous frequency for this time point, recovered with both the 
local time-frequency analysis and Hilbert transform meth- 
ods, is somewhat higher than the actual frequency. The 
fact observed here, the input data frequency and the re- 
covered instantaneous frequency being not equal, is often 
overlooked (however, see Hempelmann [2QQ2I ) . The instan- 
taneous frequency is not an entirely local value, but it also 
depends on its near neighbourhood due to the integral na- 
ture of the Hilbert transform. 

3.6. Nonharmonic components 

All the above examples were based on simple harmonic 
components. The real variable star light curves are seldom 
so simple. To show how the carrier fit method works for 
non-harmonic waveforms, we define a model with two in- 
terfering components f{t) = ci{t) -\-C2{t) -\-e{t)^ where both 
components are of the form 

c{t) = max(cos(27rti/) , 0) , (23) 

with the frequencies ui = 0.148214 and 1^2 = 0.149214, 
respectively, corresponding to the periods used in the ex- 
ample of two beating waves. Sect. 13.41 From this example, it 
is useful to recall that the beat period of these two frequen- 
cies is roughly 2000 days. In this model, however, the two 
waveforms are "active" only during a half of their cycles. 

It is quite clear that for the modelling of the process we 
need higher carrier harmonics. In Fig. [15] we show the re- 
sults of the carrier fit with K = 3. Again, a regular pattern 
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Fig. 16. Blazhko effect, phase diagram, bi-periodic fit with 
Po = 0.397680 {K = 5), P5 = 16.5462 (L = 10). See 
Sect. [3?fl 



of minima and maxima is visible. If we compare with the 
simple beat pattern of Fig. [9l the most important difference 
is the occasional bi-modality of the phase curves. This is a 
result of an interplay between phases of the two "activity 
zones" . If the zones occur at phases nearby each other, we 
have curves with a single maximum; in case of roughly 180° 
separation, two maxima appear. For more complex wave- 
forms, the interpretation of the phase plots can be even 
more convoluted. 

3. 7. Blazhko effect 

As an example of an application of the carrier fit algorithm, 
based on trigonometric modulation curves, to real observa- 
tional data, we re-analyse a subset of the photometry for 
the RR Lyr-type star MW Lyr, collected by Jurcsik et al. 
(2008). This pulsating star is known to exhibit a modula- 
tion of its light curve in addition to the variability due to 
the pulsations; this effect is known as the Blazhko effect. 
Jurcsik et al. (j2008p proposed a bi-periodic model with the 
very accurate pulsational period of Po=0. 397680 days mod- 
ulated by a significantly longer period of P5=16.5462 days. 

We start from a similar bi-periodic model with the fre- 
quencies kuo + lub^ where /c = 1, . . . , 5 and / = —10, . . . , 10, 
and uq = 1/Po and = I/P5. The main difference to the 
original analysis by Jurcsik et al. (2008) is the usage of a 
smoother form of the basic waveforms {\k\ < 5 instead of 
\k\ < 13), and a more detailed modulation model (|/| < 10 
instead of |/| < 4). The resulting normalized phase diagram 
is presented in Fig. [161 

To apply the carrier fit method to the same input data 
sequence, we fixed the "data period", D, to be 150 days 
with the coverage factor of C ^ 1.14, so that the trigono- 
metric modulation curves cover slightly more than the full 
data span. It is important to state that the resulting phase 
plot is not very sensitive to this number. We also set the 
number of harmonics = 5 for the base frequency, and 
L = 10 for the modulations. In this way the system of lin- 
ear equations to be fitted into the data was equal to the 
number of equations in the bi-periodic fit. The results are 
presented in a phase diagram of Fig. [TTl showing a strik- 
ing similarity of the phase patterns to the ones plotted in 
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Fig. 17. Blazhko effect, phase diagram, carrier fit with 
Po = 0.397680, L>=150, = 5, and L = 10. To be com- 
pared with Fig. [T6l 

Fig. [161 An important point here is that the characteristic 
modulation pattern is revealed without pre-set (or precom- 
puted) value for the modulation period. The only difference 
between the two diagrams is the somewhat larger number 
of small details for the bi-periodic fit and its exact peri- 
odicity. The phase diagram for the carrier fit solution is 
somewhat smoother and reveals minor fluctuations in the 
general periodic pattern. Whether these are fluctuations re- 
sulting form the observational noise, method artefacts, or 
sampling irregularities, remains currently open. 

As this example clearly demonstrates, it is possible to 
recover the overall bi-periodic phase structure even without 
knowing the value of one of the periods. 

3.8. Flip-flop 

To demonstrate the spline method in action, we performed 
both a trigonometric and spline-based fit into a subset of 
photometry of the star FK Com. A flip-flop phenomenon 
has earlier been detected in this very subset of data (Jetsu 
et al. '1993) covering one single observing season, during 
which a sharp jump of roughly 0.5 in phases occurs. The 
results of our analysis are shown in Figs. [T8l and [T9| for both 
the trigonometric and spline fits, the general picture is quite 
similar, and the abrupt jump in phase is clearly revealed. 
There are, however, some subtle differences. The trigono- 
metric method is in a sense more "global", due to which 
sharp features in one location can show up as oscillations 
somewhere far away. The spline method is essentially more 
"local" , and as a result the picture is somewhat smoother. 
We deliberately chose our parametrization here so that the 
total number of free parameters is equal for both methods. 

3.9. Computation of the carrier frequency 

Finally, we illustrate the method to estimate the carrier 
frequency. Our input data model is the same as used in 
Sect. 13.31 Eq. (p!9]) . containing an abrupt change in the pe- 
riod. The noise added to the input data through the term 
e{t) is now increased to 20% level. A fragment of a standard 
phase-dispersion spectrum (Stellingwerf 11978]) for this data 
set is depicted in Fig. |20l The strongest peak is at a period 




Time (days) 



Fig. 18. FK Com flip-flop, phase diagram, carrier fit with 
trigonometric modulations, carrier uq = 2.40025, K = 3^ 
L = 3, AT = 49. See Sect. ISTSl 
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Fig. 19. FK Com flip-flop, phase diagram, carrier fit with 
spline modulations, carrier uq = 2.40025, = 3, L = 5, 
N = 49. To be compared with Fig. [THl 



Po = 6.7470, coinciding closely with Pi in the input data 
model. This is to be expected, as the input data consists of 
a long fragment with this particular period. We proceed by 
using it as a carrier period value to build a time-dependent 
phase diagram depicted in Fig. [211 This carrier value pro- 
duces a nearly horizontal stripe with some oscillatory be- 
haviour for the part of the data having a period close to 
the used carrier, and a strongly tilted pattern for the part 
with a shorter period, already familiar from Sect. 13.21 for 
the case of a mismatched carrier frequency. The situation 
can be improved by refining the carrier frequency with a 
procedure described in Sect. 12.51 using the statistic D'^ de- 
fined by Eq. (p!7|) . measuring the speed of changes in the 
phase diagram. For this example calculation, we plot this 
statistics in Fig.l22l showing a well-defined minimum at the 
frequency uq = 0.14850 = 1/6.7334. Using this frequency 
as the carrier, a time-dependent phase plot with minimum 
changes can be obtained, and is depicted in Fig. I23l Now 
the decreasing and increasing fragments of the phase dia- 
gram are well balanced. 
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Fig. 20. The Stelhngwerf's statistic (criterion) for the in- 
put data with an abrupt change in the period (the same 
model as used in Sect. [3T3|) . the strongest minimum is at 
Pq = 6.7470. Note the over ah complexity of the spectrum. 
See the discussion in Sect. 13.91 
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Fig. 21. Abrupt change in the period, phase diagram, Pi = 
6.747017, P2 = 6.7018004, carrier Pq = 6.7470, K = 1, 
L = 8. See the discussion in Sect. 13.91 



Here we again note that the optimum carrier value is 
just a formal construct, and not based on any physical prin- 
ciple or insight. However, the procedure can be used as the 
last resort if no a priori information is available. 



4. Related work 

Next we discuss the proposed algorithms in the context 
of the other well-known analysis methods for variable star 
light curves. 

4.1. The 0-C method 

In the classical 0-C (observed versus computed) method 
(see for instance Sterken, I2005p , the input light curve is 
divided into short fragments, so that for each fragment it is 
possible to securely determine a certain light curve event, 
most often a local minimum or maximum. For each event, 
a difference between the observed (O) and the computed 
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Fig. 22. Phase-change spectrum, statistic (see Eq. [T7|) 
is used as the criterion to estimate stationarity of the phase 
flow. See the discussion in Sect. 13.91 
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Fig. 23. Abrupt change in the period, phase diagrams. 
Pi = 6.747017,P2 = 6.7018004, carrier Pq = 6.7334, 
K = 1,L = 8. See the discussion in Sect. 13.91 

(C) time moments is computed. The computed values are 
commonly just values from a periodic sequence 

Te=To^ PE, (24) 

where P is a period estimated beforehand, and E is an 
integer count of the full phases from the starting epoch Tq. 
The differences O — C are then plotted for a full time span 
of the input data, and a proper model for the deviations is 
then sought for. Most often polynomial models are used. 

In the carrier fit method proposed above, we are not 
very far off from this traditional scheme. The carrier period, 
^0 = ^(7^5 is just an analogue of the period P in the equa- 
tion for computed time point values. The fitting procedure 
itself tries to tightly model local fragments. In doing so, 
the local phase properties of the carrier wave indicate pe- 
culiarities of the original curve (minima and maxima among 
others). And finally, the smooth modulation curves are just 
analogues of the traditional polynomial models for O — C 
differences. In principle, it is even possible to modify the 
carrier fit method so that instead of trigonometric or spline 
models, simple polynomials are used as modulation curves. 

The O — C method is often used to refine the already 
roughly estimated periods. As shown above, a mismatch 
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in the carrier frequency reveals itself as a linear trend in 
the phase diagram. Therefore, to estimate the needed cor- 
rection to the carrier period, the linear model for carrier 
modulations is considered sufficient. 

The carrier fit method even shares the well-known prob- 
lems of the O — C method. Most importantly - the wrong 
(or not exact enough) base period P can totally spoil the 
O — C analysis. This is also true when the choice of the 
proper carrier frequency is considered. 

4.2. Fourier analysis 

Methods based on, or connected to, standard Fourier anal- 
ysis are certainly the most popular in variable star research. 
The classical series of papers, starting from Earning (1963), 
and ending with Koen (1990) and Frescura et al. (2008), 
gives a statistically sound treatment for the case when the 
input data can be described by a simple harmonic. The 
generalizations for non-harmonic (but periodic) and mul- 
tiperiodic cases are obvious (see the chapter "Frequency 
Analysis" and references in the monograph by Aerts et 
al. I2010p . For us the important point is that the classi- 
cal methods are based on the assumption that the Fourier 
spectra of the underlying processes are discrete, and form 
a sparse set of frequencies. The sparseness of the spectrum 
makes it possible to recover the constituent waveforms from 
data sets with long gaps. Unfortunately, the periodicities 
caused by the gaps can be a source of serious problems. 
Along with every discrete peak in the original spectrum 
there will be a set of spurious peaks, and a perfect recovery 
can be a challenge. 

Another method of frequency analysis is based on the 
idea of inversion. In this type of methods, the full dis- 
cretized Fourier spectrum is fitted into data (see for in- 
stance Kuhn [T9821 or more recently Nygren & Ulich [20T0| . 
For a sensible inversion, the input data needs to be sam- 
pled quite densely. It must also be noted that inversion 
techniques involve huge matrix computations, and for this 
reason these methods fail to work for long input data sets. 

The carrier fit algorithm described above combines 
properties of the discrete and continuous spectrum meth- 
ods. First, we have a carrier frequency and its overtones 
kvQ^ which constitutes the discrete aspect. If the modu- 
lating waveforms ak(t) and hkit) are smooth enough, the 
resulting spectrum will consist of K narrow but continu- 
ous frequency bands, accounting for the continuous aspect. 
By combining some good properties of the standard meth- 
ods, we can achieve a compromise, where the discreteness 
of the midpoints of the narrow bands allow to "bridge" long 
gaps in time domain, and the continuity of the local spec- 
tra around the mid-frequencies allows to describe smooth 
time-dependent changes. Our method can essentially be re- 
garded as a quite straightforward combination of the two 
classical approaches. 

4.3. Local frequency analysis 

There is another set of methods popular among variable 
star researchers. These methods are based on a division of 
the input data into separate segments, and then using some 
standard method to analyse them. The segmentation can 
be based on natural fragmentation of the input data or by 
formal division of the full data span. In some methods the 



segments can be overlapping in time (see recent detailed 
implementation by Lehtinen et al 2011). It is assumed that 
the changes in the local Fourier spectra or other type statis- 
tics computed for different segments reveal actual trends in 
frequency structure of the underlying signal. 

A typical example - the 5 Scuti -type star Tuc was 
observed during six nights (Stobie & Shobbrook 1976), and 
data for each night was Fourier analysed seprately. The 
conclusion of the authors was 

Photometric observations of Tucanae were Fourier 
analysed for their component frequencies. No coher- 
ent frequency could be found spanning the complete 
data set. It is shown that both the frequencies and 
amplitudes present in Tucanae change on time 
scale as short as 24 hr. 

However, more recent analysis of more extended data sets 
by other authors (for example Kurtz I1980[ Paparo et 
al. I1996P gave another interpretation to the local variability 
of Tuc. It occurred that the varying local spectra result 
from constructive and destructive interferences of different 
coherent oscillation modes. The interpretation of the lo- 
cal spectra is further complicated by their low resolution. 
Above we gave an illustration of this effect (see Fig. [TQ|) 
in the context of straightforward local Fourier analysis. 
This example shows that it is a good practice to use time- 
dependent phase diagrams (built by the carrier fit method) 
and time-dependent frequency analysis together to get a 
full understanding of the underlying process. 

Localisation is carried to the limit in the method pro- 
posed by Tsantilas & Rovithis-Livaniou (|2008p . They use a 
very similar light curve model to Eq. ([T]) 

f(t) = a{t) sm{b{t) X t + c(t)), (25) 

where observed data is decomposed into three continu- 
ous curves: amplitude a(t), frequency b{t) and phase c{t). 
The point estimates for the three components are then ob- 
tained by performing local non-linear least squares fits. 
Unfortunately, the presentation of the light curve in this 
form is not unique (for instance frequency changes can be 
compensated by phase changes) and consequently the ob- 
tained numerical results contain large amount of arbitrari- 
ness. In our method the carrier is fixed and modulations are 
constrained to be smooth. This allows to get more useful 
presentations for the light curves. 

4.4. Time-frequency analysis 

Another set of methods worthwhile to mention in this con- 
text are different time- frequency transforms. They are all 
based on a decomposition of the input data into set of sub- 
waves. Every sub-wave (in some methods called wavelet) 
is localized in time as well as in frequency. For a large set 
of this kind of decompositions we refer to a recent work 
by Blackman (|201Q| and references therein. In the context 
of variable star research, the most important problem with 
the standard time-frequency methods is related to the gaps 
in the input data. Usually the set of sub-waves is mathe- 
matically complete, but this leads to unrealistic "filling of 
gaps" . The sub- waves which are localized in short time in- 
tervals inside observing gaps will get zero values and as a 
result of that the structure of gaps will be seen through 



12 



J. Pelt et al: Carrier fit method 



in the final two-dimensional spectrum. A very nice illustra- 
tion of this point can be found in Szatmary (1994), where 
a large catalogue of model wavelet decompositions is given. 

To overcome problems with gaps it is possible to modify 
the original wavelet method by making its resolution to 
depend on time. This can be achieved by the use of certain 
weighting schemes (see Foster 1996 for details) or adaptive 
wavelets, as done by Frick et al. (1997). 

In the carrier fit method we use a priori information to 
constrain our sub-waves to be localised in fixed frequency 
strips around carrier and possibly around its harmonics. 
The shortcoming of our formulation is that our system 
of sub- waves is not complete. Nevertheless, with our al- 
gorithm, the gaps can be properly filled in, which enables 
to get a continuous phase or frequency diagrams that are 
more easy to interpret. 

5. Discussion and conclusions 

5.1. Utility of the carrier fit 

All the simple examples above demonstrate that the car- 
rier fit procedure itself is not sufficient to perform proper 
data analysis and interpretation. It is just a particular 
method to interpolate data sets with gaps to obtain a rel- 
atively smooth estimate for the variable star light curve. 
Interpolation itself has some value only when the real ob- 
served process and its observed values have certain prop- 
erties. Most importantly, the process must be governed 
by some relatively stationary inherent clocking mechanism 
whose parameters change slowly. Fortunately this is the 
case for many variable stars. 

5.2. Domain of applicability 

By formulation and by algorithmic implementation the pro- 
posed utility can be considered as automatic. However, the 
obtained results can be very often misleading. This is why 
the complete understanding of the involved ideas is very 
important. Let us now reiterate major requirements for suc- 
cessful application of the carrier fit method. 

Most importantly, the data to be analysed must come 
from an object with a certain internal clocking mechanism, 
being it rotation, orbiting or pulsation. This allows to fix, at 
least approximately, a proper value for a carrier frequency. 

The modulating curves of the basic cycle must be rela- 
tively slow so that the spectral representations of the modu- 
lations and carrier frequency are located in different parts of 
the spectrum, i.e. the Bedrosian theorem is valid. It is also 
important that models for modulations must be smooth 
enough to allow "bridging the gaps" in observational se- 
quences. 

There must be enough of high quality observational 
points to allow adequate interpolation. To avoid over fit 
and under fit, standard methods from regression analysis 
should be used (Draper & Smith I1998|) . 

The user of the methods must have a certain expertise 
to interpret the obtained phase diagrams. The proper and 
required intuition for that can be developed by using model 
computations with known input components. As a starting 
point the examples provided above can be useful. 

Because the method involves inversion of large matri- 
ces, certain numerical problems may occur. The methods 
based on row-wise updating (Givens rotations) to solve the 



least-squares computational problems are very useful in this 
context; for details see Lawson & Hanson (I1987p . 

The carrier fit method belongs to the class of data 
analysis methods called explorative. After getting an in- 
tuitively pleasing and convincing phase diagrams the re- 
searcher must quantify the results using more precise mod- 
els (with proper estimates of errors etc). 

Another important aspect of the new method is the fo- 
cus on phase behaviour. For that purpose, the visualization 
of normalized light curves is of great importance. This facil- 
itates the detection of relevant trends and abrupt changes 
in light curves. This is just the aspect which usually re- 
mains hidden when standard multicomponent fit into the 
data is used. 

In this paper we have presented a relatively new but 
simple class of data analysis methods which are useful in 
different astrophysical contexts. The best way to promote 
the methods would be their successful application to actual 
data sets. We are undertaking this effort in the next papers 
of the series. 

The developed software is a part of the general time 
series package (Irregularly Spaced Data Analysis) which is 
freely available q 
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